Analysis of Elastodynamic Deformations near a Crack/Notch Tip by the Meshless Local Petrov-Galerkin (MLPG) Method

The Meshless Local Petrov-Galerkin (MLPG) method is used to analyze transient deformations near either a crack or a notch tip in a linear elastic plate. The local weak formulation of equations governing elastodynamic deformations is derived. It results in a system of coupled ordinary differential equations which are integrated with respect to time by a Newmark family of methods. Essential boundary conditions are imposed by the penalty method. The accuracy of the MLPG solution is established by comparing computed results for one-dimensional wave propagation in a rod with the analytical solution of the problem. Results are then computed for the following two problems: a rectangular plate with a central crack with plate edges parallel to the crack axis loaded in tension, and a double edge-notched plate with the edge between the notches loaded by compressive tractions. Stresses at points near the crack/notch tip computed from the MLPG solution are found to agree well with those obtained from either the analytical or the finite element solution of the same problem. The index of stress singularity is ascertained from a plot of log (stress) vs. log ( r) where r is the distance from the crack tip. It is found that, for the double-edge notched plate, the mode-mixity of deformations near a notch-tip in an orthotropic plate can be adjusted by suitably varying the in-plane moduli of the material of the plate. The variation of shear stress with r exhibits a boundary layer effect near r= O. tegrals appearing in the local weak formulation of the problem. Atluri et al. (1999) have pointed out that the Galerkin approximation can also be adopted that leads to a symmetric stiffness matrix. Atluri and Zhu (2000) solved elastostatic problems by the MLPG method, and Lin and Atluri (2000) introduced the upwinding scheme to analyze steady convection-diffusion problems. Ching and Batra (200 1) enriched the polynomial basis functions with those appropriate to describe singular deformation fields near a crack tip and used the diffraction criterion to find stress intensity factors, the J-integrals and singular stress fields near a crack tip. Gu and Liu (2001) used the Newmark family of methods to study forced vibrations of a beam. The problem of bending of a thin plate has been studied by Long and Atluri (2002). Warlock et al. (2002) have analyzed elastostatic deformations of a material compressed in a rough rectangular cavity analytically by the Laplace transformation technique and numerically by the MLPG method. Atluri and Shen (2002a,b) have demonstrated the use of different weight functions and have compared their performance with that of the Galerkin finite element method. By choosing a Heaviside step function as the test function, they eliminated the domain integration in the local weak form. Thus only boundary integrals over local subdomains remained in the local weak form. For elastostatic problems, this was shown to be more efficient than both the finite element and the boundary element methods. The paper is organized as follows. Section 2 gives the MLPG formulation including the local weak form, the moving least squares approximation, the discrete governing equations and the time integration scheme. Calculations of the dynamic stress intensity factors from the near-tip stress fields are also described. Numerical examples are presented in Section 3. The MLPG results are compared with either analytical or finite element solutions. Section 4 summarizes the conclusions. Introduction The meshless method has attracted considerable attention in the past two decades due to the flexibility of placing nodes in the domain of study. Atluri and Zhu (1998) have proposed a meshless method which requires no background mesh to evaluate numerically various inI Department of Engineering Science and Mechanics, MC 0219 Virginia Polytechnic Institute and State University RI"...J,.~hI,ra VA '),101;1 718 Copyright @ 2002 Tech Science Press 2 Formulation of the Problem Integrating the first term on the left side of (8) by parts, and using natural boundary condition (7) we obtain 2.1 Governing equations f pv iUid.Q. + f v i,jO"ijd.Q. + a f J.o.s J.o.s Jrsu VjUjdr For a plane linear elastodynarnic problem on domain .Q. bounded by the boundary r, governing equations in rectangular Cartesian coordinates are Vitidr = r vJidr+a r JrS! Jrsu ViUidr + { vibid.Q, lo., (9) O'ij,j+bi=pui,in.o.,(i,j=1,2), (1) O'ij = ~kkOij + 2,uEij, in.o., (2) Eij=(ui,j+uj,i)/2,in.o.. (3) Here P is the mass density, Ui the displacement, t the time, Ui = azui/dt2 the acceleration, O'ij the stress tensor, Eij the infinitesimal strain tensor, b i density of the body force vector, A and,u are Lame' constants for the material of the body, Ui,j = dUi/dXj, x gives the present position of a material particle, and a repeated index implies summation over the range of the index. Equations (1)-(3) are supplemented with the following initial and boundary conditions: (4) where the unknown coefficients a(x,t) are functions of the space coordinates xT = [XI,X2] and time t, and p(x) is (5) a complete monomial in x having m terms. The complete (6) quadratic monomials basis functions in two-dimensions (7) are u(x,to) = UQ(x), x E .0., ti(X,tO) = tio(X), x E .0., Ui = Ui, on r u, ti == O'ijnj = ii, on rto Here Uj,tj,no and Do denote the prescribed displacepT(x)=[1,Xl,X2,(Xl)2,XlX2,(x2)2j;m=6. (11) ments, tractions, initial displacements and initial velocities, respectively, n j is the unit outward normal to r, and For each component of u, the coefficients a(x,t) in (10) r u and r t are complementary parts of r where essential are obtained by minimizing J defined by and natural boundary conditions are prescribed. n J = LW(X-Xi)[pT (Xi)a(X,t) Ui(t)]2 ---' Implementation of the MLPG method i=l Taking the inner product of (I) with v and of (6) with av, Here Ui is the fictitious value of a component of u h at and integrating the resulting equations over .0. s and r su x = Xi, and n is the number of nodes in the domain of respectively, we obtain influence of x for which the weight functions w( x x J # O. Several different weight functions are given in Atluri 1 (O'ij,j pili + bi)Vid.o. -1 a(Ui Ui)Vidr = 0, (8) and Shen (2002a,b); here, the following Gaussian weight Qs r su function is used. 2.2 where r su = r s n r u, r s is the boundary of the local domain .o.s c .0. and a is a penalty parameter used to satisfy the essential boundary conditions. For eqn. (8) to be dimensionally correct, a must have units of Force/(Length)3. The penalty parameter a may vary from point to point but is usually taken to be a constant with magnitude much larger than 1../ L where L is a typical dimension of the body. Henceforth we take a to be a constant. [ )2k] 2k r' J!=!i.l exp ~ exp [ ( Cj ) ] (Cj , 0:$ Ix xii :$ Ti l-eXp[=(~)2k] W(X-Xi) =